# This script creates all of the tables and figures in the paper

rm(list = ls())

library(tidyverse)
library(lfe)
library(estimatr)
library(stargazer)
library(patchwork)

############## Specify file paths ##############

# path to final dataset
datapath <- "[FILEPATH HERE]"

# output path
setwd("[PATH TO OUTPUT FOLDER HERE]")


############## Final Regression Setup ##############

# read final dataset
main <- readRDS(paste0(datapath, "final_data.rds"))

# set up event time as a factor
main <- main %>% 
  mutate(eventtime = -1*(fire == 0) + (FISCAL_YEAR - eventyear)*(fire == 1))
main$eventtime <- relevel(factor(main$eventtime), ref = 5)
main <- main %>% 
  mutate(treatDID = 1*(as.character(eventtime) %in% c("0", "1", "2", "3", "4")))

# inverse hyperbolic sine transformation
ihs <- function(x) {
  y <- log(x + sqrt(x ^ 2 + 1))
  return(y)
}

# add the omitted -1 year
omit <- data.frame(term = "-1", 
                   estimate = 0,
                   std.error = 0,
                   statistic = 0,
                   p.value = 0,
                   conf.low = 0,
                   conf.high = 0)


############## Figure 1 ##############

felm(ihs(POP) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
     data = main) %>% summary()

coef <- felm(ihs(POP) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term)) 

coef <- rbind(coef[1:4,], omit, coef[5:9,])

annotation <- data.frame(
  x = 2,
  y = -0.06,
  label = "DD coefficient = -0.00759 (0.012)"
)

ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  geom_segment(aes(x = -1, y = -0.0076, xend = 4, yend = -0.0076), 
               col = "red", size = 0.4, linetype = "dashed") +
  geom_text(data=annotation, aes( x=x, y=y, label=label)) +
  annotate("segment", x = 0, xend = -0.5, y = -0.05, yend = -0.01, width = 0.1, 
           arrow=arrow(length = unit(0.1, "inches"))) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) +
  ggtitle("Population (asinh)") + ylim(-0.1, 0.1) +
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))

ggsave("figure1.pdf", device = cairo_pdf, width = 6, height = 4)


############## Figure 2 ##############

# general revenue 

coef <- felm(ihs(GENREV_TOT) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term)) 

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p1 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") + 
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.2, 0.4) +
  ggtitle("A. Total (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))

# tax revenue
coef <- felm(ihs(GENREV_TOT_TAXES) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term)) 

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p2 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.2, 0.4) +
  ggtitle("B. Taxes (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))


# property tax 

coef <- felm(ihs(GENREV_PROPTAX) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term)) 

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p3 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.3, 0.6) +
  ggtitle("C. Property tax (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))


# sales tax
coef <- felm(ihs(GENREV_SALE_USE_TAX) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term)) 

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p4 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.3, 0.6) +
  ggtitle("D. Sales and use tax (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))

(p1 + p2) / (p3 + p4)

ggsave("figure2.pdf", device = cairo_pdf, width = 8, height = 6)



############## Figure 3 ##############

# total functional revenue

coef <- felm(ihs(FUNC_TOT) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term)) 

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p1 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.4, 0.8) +
  ggtitle("A. Total (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))


# total functional revenue from taxes

coef <- felm(ihs(FUNC_TOT_TAXES) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term)) 

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p2 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-3,6) +
  ggtitle("B. Special taxes (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))


# functional revenue from charges

coef <- felm(ihs(FUNC_CHARGES) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term)) 

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p3 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.5, 1.0) +
  ggtitle("C. Charges (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))


# intergovernmental transfers - functional

main <- main %>% rowwise() %>% mutate(FUNC_INT = sum(FUNC_FED, FUNC_ST, FUNC_CO, na.rm = T))

coef <- felm(ihs(FUNC_INT) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term)) 

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p4 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.5, 1.0) +
  ggtitle("D. Intergov transfers (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))

(p1 + p2) / (p3 + p4)

ggsave("figure3.pdf", device = cairo_pdf, width = 8, height = 6)


############## Figure 4 ##############

# total expenditure

coef <- felm(ihs(TOT_EXP_TOT) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p1 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.2, 0.4) +
  ggtitle("A. Total (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))

# public safety
coef <- felm(ihs(TOT_EXP_PUB_SAFETY) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p2 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.2, 0.4) +
  ggtitle("B. Public safety (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))


# general government
coef <- felm(ihs(TOT_EXP_GEN_GOV) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p3 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.5, 1.0) +
  ggtitle("C. General government (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))


# community development

coef <- felm(ihs(TOT_EXP_COMM_DEVELOP) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main %>% filter(TOT_EXP_COMM_DEVELOP > 0)) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p4 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.5, 1.0) +
  ggtitle("D. Community development (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))


# transportation
coef <- felm(ihs(TOT_EXP_TRANSP) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p5 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.8, 1.2) +
  ggtitle("E. Transportation (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))


# culture and leisure

coef <- felm(ihs(TOT_EXP_CULT_LEISURE) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p6 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.8, 1.2) +
  ggtitle("F. Culture and leisure (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))

# health
coef <- felm(ihs(TOT_EXP_HEALTH) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p7 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-2, 3) +
  ggtitle("G. Health (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))

# public utilties
coef <- felm(ihs(TOT_EXP_PUB_UTILITIES) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p8 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.8, 1.2) +
  ggtitle("H. Public utilities (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))

(p1 + p2) / (p3 + p4) / (p5 + p6) / (p7 + p8)

ggsave("figure4.pdf", device = cairo_pdf, width = 8, height = 9)



############## Figure 5 ##############

coef <- felm(identity(EXP_REV_TOT/POP) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main %>% mutate(EXP_REV_TOT = - EXP_REV_TOT)) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p1 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), 
              fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-500, 200) +
  ggtitle("A. Excess functional revenue per capita") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))

coef <- felm(identity(GENREV_EXC_OF_GEN_REV/POP) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

p2 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), 
              fill = "blue", alpha = 0.2) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-500, 200) +
  ggtitle("B. Excess total revenue per capita") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))

p1 + p2 

ggsave("figure5.pdf", device = cairo_pdf, width = 8, height = 3)



############## Figure 6 ##############


coef <- felm(DEFICIT ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main %>% mutate(DEFICIT = GENREV_EXC_OF_GEN_REV < 0)) %>% 
  tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

annotation <- data.frame(
  x = 2,
  y = -0.2,
  label = "DD coefficient = 0.247***(0.050)"
)

ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), 
              fill = "blue", alpha = 0.2) +
  geom_segment(aes(x = -1, y = 0.247, xend = 4, yend = 0.247), 
               col = "red", size = 0.4, linetype = "dashed") +
  geom_text(data=annotation, aes( x=x, y=y, label=label)) +
  annotate("segment", x = 2, xend = 2, y = -0.1, yend = 0.2, width = 0.1, 
           arrow=arrow(length = unit(0.1, "inches"))) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.5, 1.0) +
  ggtitle("1(Overall budget deficit)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))


ggsave("figure6.pdf", device = cairo_pdf, width = 6, height = 4)


############## Figure 8 ##############

coef <- felm(ihs(GENREV_REAL_PROP_TRANSFER_TAX) ~ 
               eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% 
  tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term)) 

coef <- rbind(coef[1:4,], omit, coef[5:9,])

annotation <- data.frame(
  x = 1,
  y = 1.2,
  label = "DD coefficient = 0.486* (0.265)"
)


ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  geom_segment(aes(x = -1, y = 0.486, xend = 4, yend = 0.486), 
               col = "red", size = 0.4, linetype = "dashed") +
  geom_text(data=annotation, aes( x=x, y=y, label=label)) +
  annotate("segment", x = 1, xend = 1, y = 1.1, yend = 0.5, width = 0.1, 
           arrow=arrow(length = unit(0.1, "inches"))) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-2, 2) +
  ggtitle("Property transfer tax (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))

ggsave("figure8.pdf", device = cairo_pdf, width = 6, height = 4)


############## Figure 9 ##############

# asinh total fire expenditure

coef <- felm(ihs(EXP_FIRE_PUB) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

annotation <- data.frame(
  x = 2.5,
  y = -3,
  label = "1.494*** (0.533)"
)

p1 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), 
              fill = "blue", alpha = 0.2) +
  geom_segment(aes(x = -1, y = 1.494, xend = 4, yend = 1.494), 
               col = "red", size = 0.4, linetype = "dashed") +
  geom_text(data=annotation, aes( x=x, y=y, label=label)) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-4, 6) +
  ggtitle("A. Fire (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))

# total fire expenditure level

coef <- felm(identity(EXP_FIRE_PUB/1000) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

annotation <- data.frame(
  x = 2.5,
  y = -1200,
  label = "685.7* (351.4)"
)

p2 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), 
              fill = "blue", alpha = 0.2) +
  geom_segment(aes(x = -1, y = 685.7, xend = 4, yend = 685.7), 
               col = "red", size = 0.4, linetype = "dashed") +
  geom_text(data=annotation, aes( x=x, y=y, label=label)) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-1600, 2400) +
  ggtitle("B. Fire (level in 1,000s)") +
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))


# asinh total disaster expenditure

main <- main %>% mutate(OPEXP_DISASTER_PUB = as.numeric(OPEXP_DISASTER_PUB))

coef <- felm(ihs(OPEXP_DISASTER_PUB) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

annotation <- data.frame(
  x = 2.5,
  y = -4,
  label = "3.37*** (1.593)"
)

p3 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), 
              fill = "blue", alpha = 0.2) +
  geom_segment(aes(x = -1, y = 3.37, xend = 4, yend = 3.37), 
               col = "red", size = 0.4, linetype = "dashed") +
  geom_text(data=annotation, aes( x=x, y=y, label=label)) +
  xlab("Year relative to fire") + #ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-6, 12) +
  ggtitle("C. Disaster prep. (asinh)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))

# total disaster expenditure

coef <- felm(identity(OPEXP_DISASTER_PUB/1000) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

annotation <- data.frame(
  x = 2.5,
  y = -67,
  label = "57.01*** (19.46)"
)

p4 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), 
              fill = "blue", alpha = 0.2) +
  geom_segment(aes(x = -1, y = 57, xend = 4, yend = 57), 
               col = "red", size = 0.4, linetype = "dashed") +
  geom_text(data=annotation, aes( x=x, y=y, label=label)) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + #ylim(-100, 200) +
  ggtitle("D. Disaster prep. (level in 1,000s)") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))


(p1 + p2)/(p3 + p4) + plot_annotation(theme = theme(plot.title = element_text(size = 16)))

ggsave("figure9.pdf", device = cairo_pdf, width = 8, height = 6)


############## Figure A3 ##############

# extensive margin

coef <- felm(HEALTH_IND ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main %>% mutate(HEALTH_IND = TOT_EXP_HEALTH > 0)) %>% 
  tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

annotation <- data.frame(
  x = 1.5,
  y = -0.05,
  label = "0.0786*** (0.019)"
)

p1 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  geom_segment(aes(x = -1, y = 0.0786, xend = 4, yend = 0.0786), 
               col = "red", size = 0.4, linetype = "dashed") +
  geom_text(data=annotation, aes( x=x, y=y, label=label)) +
  annotate("segment", x = 1.5, xend = 1.5, y = -0.03, yend = 0.06, width = 0.1, 
           arrow=arrow(length = unit(0.1, "inches"))) +
  xlab("Year relative to fire") + ylab("Estimate") + 
  scale_x_continuous(breaks = c(-5:4)) + ylim(-0.2, 0.2) +
  ggtitle("A. Extensive margin") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))

# intensive margin

coef <- felm(ihs(TOT_EXP_HEALTH) ~ eventtime | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main %>% filter(TOT_EXP_HEALTH > 0)) %>% tidy(conf.int = TRUE)

coef$term <- as.character(coef$term) %>% substr(10, nchar(coef$term))

coef <- rbind(coef[1:4,], omit, coef[5:9,])

annotation <- data.frame(
  x = 1.5,
  y = -0.8,
  label = "-0.254 (0.187)",
  size = 4)

p2 <- ggplot(coef, aes(x = as.numeric(term), y = estimate)) + 
  geom_hline(yintercept = 0, col = "gray") + 
  geom_vline(xintercept = -1, linetype = "dashed", alpha = 0.7) +
  geom_point(aes(y = estimate), size = 1, color = "blue") + 
  geom_line(aes(y = estimate), color = "blue") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue", alpha = 0.2) +
  geom_segment(aes(x = -1, y = -0.254, xend = 4, yend = -0.254), 
               col = "red", size = 0.4, linetype = "dashed") +
  geom_text(data=annotation, aes( x=x, y=y, label=label)) +
  annotate("segment", x = 1.5, xend = 1.5, y = -0.7, yend = -0.3, width = 0.1, 
           arrow=arrow(length = unit(0.1, "inches"))) +
  xlab("Year relative to fire") + ylab("Estimate") +
  scale_x_continuous(breaks = c(-5:4)) + ylim(-1.5, 1) +
  ggtitle("B. Intensive margin") + 
  theme_classic() +
  theme(panel.grid.major.y = 
          element_line(size = 0.5, linetype = 'solid', colour = "gray"))


p1 + p2 + plot_annotation(theme = theme(plot.title = element_text(size = 14)))

ggsave("figureA3.pdf", device = cairo_pdf, width = 8, height = 3)


############## Table 1 ##############

# read municipality-year panel
munipanel <- readRDS(paste0(datapath, "muni_panel.rds"))

munipanel <- munipanel %>% rowwise() %>%
  mutate(FUNC_INT = sum(FUNC_FED, FUNC_ST, FUNC_CO, na.rm = T)) %>%
  mutate(EXP_REV_TOT = -EXP_REV_TOT,
         A = (SERVICE == "a")*100,
         B = (SERVICE == "b")*100,
         C = (SERVICE == "c")*100,
         D = (SERVICE == "d")*100,
         E = (SERVICE == "e")*100,
         F = (SERVICE == "f")*100,
         X = (SERVICE == "x")*100,
         CHARTER = (CLASS_DATA == "Chartered")*100,
         DEFICIT = (GENREV_EXC_OF_GEN_REV < 0)*100)


# list of treated cities
treated <- munipanel %>% 
  filter(percent_fire > 0.10, FISCAL_YEAR >= 1995) %>% 
  select(ENTITY_ID, FISCAL_YEAR)


treated_p <- munipanel %>%
  filter(ENTITY_ID %in% treated$ENTITY_ID) %>% 
  select(POP, GENREV_TOT, GENREV_TOT_TAXES, GENREV_PROPTAX, GENREV_SALE_USE_TAX,
         GENREV_REAL_PROP_TRANSFER_TAX, FUNC_TOT, FUNC_TOT_TAXES, FUNC_CHARGES, 
         FUNC_INT, TOT_EXP_TOT, TOT_EXP_PUB_SAFETY, TOT_EXP_GEN_GOV,  
         TOT_EXP_COMM_DEVELOP, TOT_EXP_TRANSP, TOT_EXP_CULT_LEISURE, 
         TOT_EXP_HEALTH, TOT_EXP_PUB_UTILITIES, EXP_FIRE_PUB, EXP_DISASTER_PUB,
         EXP_REV_TOT, GENREV_EXC_OF_GEN_REV, DEFICIT, CHARTER,
         A, B, C, D, E, F, X) %>%
  data.frame()

for (i in 2:22) {
  treated_p[, i] <- treated_p[, i]/treated_p$POP
}

stargazer(treated_p, type = "text", out = "Table1_part1.txt", 
          summary.stat = c("n", "mean", "sd"), digits = 1)

untreated_p <- munipanel %>% rowwise() %>%
  mutate(FUNC_INT = sum(FUNC_FED, FUNC_ST, FUNC_CO, na.rm = T)) %>%
  mutate(EXP_REV_TOT = -EXP_REV_TOT,
         A = (SERVICE == "a")*100,
         B = (SERVICE == "b")*100,
         C = (SERVICE == "c")*100,
         D = (SERVICE == "d")*100,
         E = (SERVICE == "e")*100,
         F = (SERVICE == "f")*100,
         X = (SERVICE == "x")*100,
         CHARTER = (CLASS_DATA == "Chartered")*100,
         DEFICIT = (GENREV_EXC_OF_GEN_REV < 0)*100)

untreated_p <- untreated_p %>%
  filter(!ENTITY_ID %in% treated$ENTITY_ID,
         ENTITY_NAME != "San Francisco") %>% 
  select(POP, GENREV_TOT, GENREV_TOT_TAXES, GENREV_PROPTAX, GENREV_SALE_USE_TAX,
         GENREV_REAL_PROP_TRANSFER_TAX, FUNC_TOT, FUNC_TOT_TAXES, FUNC_CHARGES, 
         FUNC_INT, TOT_EXP_TOT, TOT_EXP_PUB_SAFETY, TOT_EXP_GEN_GOV,  
         TOT_EXP_COMM_DEVELOP, TOT_EXP_TRANSP, TOT_EXP_CULT_LEISURE, 
         TOT_EXP_HEALTH, TOT_EXP_PUB_UTILITIES, EXP_FIRE_PUB, EXP_DISASTER_PUB,
         EXP_REV_TOT, GENREV_EXC_OF_GEN_REV, DEFICIT, 
         CHARTER, A, B, C, D, E, F, X) %>%
  data.frame()

untreated_p$POP[untreated_p$POP==0] <- NA

for (i in 2:22) {
  untreated_p[, i] <- untreated_p[, i]/untreated_p$POP
}

stargazer(untreated_p, type = "text", out = "Table1_part2.txt", 
          summary.stat = c("n", "mean", "sd"), digits = 1)



############## Table 2 ##############

main1 <- main %>% filter(!is.na(obs)) %>%
  mutate(white = white*100,
         black = black*100,
         hisp = hisp*100,
         occup = occup*100)


reg1 <- felm(ihs(income) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | ENTITY_ID, 
             data = main1, weights = main1$obs) 

reg2 <- felm(white ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | ENTITY_ID, 
             data = main1, weights = main1$obs) 

reg3 <- felm(hisp ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | ENTITY_ID, 
             data = main1, weights = main1$obs)

reg4 <- felm(black ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | ENTITY_ID, 
             data = main1, weights = main1$obs)

reg5 <- felm(occup ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | ENTITY_ID, 
             data = main1, weights = main1$obs) 


stargazer(reg1, reg2, reg3, reg4, reg5, type = "text", out = "table2.txt",
          title = "DD Estimates on Demographic Changes",
          dep.var.labels = c("log(Income)", "%White", "%Hispanic", "%Black", "%Occup"))

rm(main1)


############## Table 3 ##############

did1 <- felm(ihs(GENREV_TOT) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main)

did2 <- felm(ihs(GENREV_TOT_TAXES) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main)

did3 <- felm(ihs(GENREV_PROPTAX) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main)

did4 <- felm(ihs(GENREV_SALE_USE_TAX) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main)

stargazer(did1, did2, did3, did4, type = "text", out = "table3.txt",
          title = "DD Estimates of General Revenue Categories",
          dep.var.labels = c("log(GenRev)", "log(Taxes)", "log(PropTax)", "log(SalesTax)"))


############## Table 4 ##############

did1 <- felm(ihs(FUNC_TOT) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) 

did2 <- felm(ihs(FUNC_TOT_TAXES) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) 

did3 <- felm(ihs(FUNC_CHARGES) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) 

did4 <- felm(ihs(FUNC_INT) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) 

stargazer(did1, did2, did3, did4, type = "text", out = "table4.txt",
          title = "DD Estimates of Functional Revenue Categories",
          dep.var.labels = c("log(FuncTot)", "log(Taxes)", "log(Charges)", "log(InterGov)"))


############## Table 5 ##############

did1 <- felm(ihs(TOT_EXP_TOT) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) 

did2 <- felm(ihs(TOT_EXP_PUB_SAFETY) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) 

did3 <- felm(ihs(TOT_EXP_GEN_GOV) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) 

did4 <- felm(ihs(TOT_EXP_COMM_DEVELOP) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main %>% filter(TOT_EXP_COMM_DEVELOP > 0)) 

did5 <- felm(ihs(TOT_EXP_TRANSP) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) 

did6 <- felm(ihs(TOT_EXP_CULT_LEISURE) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) 

did7 <- felm(ihs(TOT_EXP_HEALTH) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) 

did8 <- felm(ihs(TOT_EXP_PUB_UTILITIES) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) 

stargazer(did1, did2, did3, did4, did5, did6, did7, did8,
          title = "DD Estimates of Expenditures", type = "text", out = "table5.txt",
          dep.var.labels = c("Total", "Public Safety", "General Government", "Community Development",
                             "Transportation", "Culture and Leisure", "Health", "Public Utilities"))


############## Table 6 ##############

did1 <- felm(ihs(EXP_FIRE_PUB) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event,
             data = main)

did2 <- felm(identity(EXP_FIRE_PUB/1000) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event,
             data = main)

did3 <- felm(ihs(OPEXP_DISASTER_PUB) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event,
             data = main)

did4 <- felm(identity(OPEXP_DISASTER_PUB/1000) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event,
             data = main)

stargazer(did1, did2, did3, did4, type = "text", out = "table6.txt",
          dep.var.labels.include = FALSE,
          column.labels   = c("Fire", "Disaster Preparedness"),
          column.separate = c(2, 2),
          title = "DD Estimates of Expenditures on Fire",
          add.lines = list(c("%Change", "286", "75.5", "721", "164"),
                           c("D.V. mean", "6.74", "908.6", "4.45", "34.77"),
                           c("Municipal FE", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE", "Yes", "Yes", "Yes", "Yes")),
          omit.stat = c("ser", "adj.rsq"))


############## Table 7 ##############

did1 <- felm(identity(EXP_REV_TOT/POP) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main %>% mutate(EXP_REV_TOT = - EXP_REV_TOT)) 

did2 <- felm(identity(GENREV_EXC_OF_GEN_REV/POP) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
             data = main) 

stargazer(did1, did2, type = "text", out = "table7.txt",
          title = "DD Estimates of Expenditures on Budget Balance",
          dep.var.labels = c("Exc. Func. Rev.", "Exc. Tot. Rev."))


############## Table 8 ##############

pop <- felm(ihs(POP) ~ treatDID | ENTITY_ID + FISCAL_YEAR + event | 0 | event, 
            data = main)

genrev <- felm(ihs(GENREV_TOT) ~ treatDID | ENTITY_ID + FISCAL_YEAR + event | 0 | event, 
               data = main)

funcrev <- felm(ihs(FUNC_TOT) ~ treatDID | ENTITY_ID + FISCAL_YEAR + event | 0 | event, 
                data = main) 

exp <- felm(ihs(TOT_EXP_TOT) ~ treatDID | ENTITY_ID + FISCAL_YEAR + event | 0 | event, 
            data = main) 

excrev <- felm(identity(GENREV_EXC_OF_GEN_REV/POP) ~ treatDID | ENTITY_ID + FISCAL_YEAR + event | 0 | event, 
               data = main) 

stargazer(pop, genrev, funcrev, exp, excrev, 
          type = "text", out = "table8.txt",
          title = "DD Estimates of Major Categories (Event FE)",
          dep.var.labels = c("Population", "Gen. Rev.", "Func. Rev.", "Total Exp.", "Tot. Exc. Rev."))



############## Table 9 ##############

weights <- main %>% 
  filter(fire == 1, as.character(eventtime) == "0") %>%
  select(event, TREATPOP = POP)

main1 <- main %>% left_join(weights, by = "event")

pop <- felm(ihs(POP) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
            data = main1, weights = main1$TREATPOP)

genrev <- felm(ihs(GENREV_TOT) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
               data = main1, weights = main1$TREATPOP)

funcrev <- felm(ihs(FUNC_TOT) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
                data = main1, weights = main1$TREATPOP)

exp <- felm(ihs(TOT_EXP_TOT) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, 
            data = main1, weights = main1$TREATPOP)

excrev <- felm(identity(GENREV_EXC_OF_GEN_REV/POP) ~ treatDID | ENTITY_ID + FISCAL_YEAR + event | 0 | event, 
               data = main1, weights = main1$TREATPOP)


stargazer(pop, genrev, funcrev, exp, excrev, type = "text", out = "table9.txt",
          title = "DD Estimates from Population-Weighted Regressions",
          dep.var.labels = c("Population", "Gen. Rev.", "Func. Rev.", "Total Exp.", "Tot. Exc. Rev."))

rm(weights, main1)



############## Table 10 ##############

# read municipality-year panel
munipanel <- readRDS(paste0(datapath, "muni_panel.rds"))

# treated dataset
treated <- munipanel %>% 
  filter(popfire > 5000,
         FISCAL_YEAR >= 1995) %>% 
  arrange(desc(percent_fire)) %>%
  filter(ENTITY_NAME != "Los Angeles" | as.character(FISCAL_YEAR) == "2008",
         ENTITY_NAME != "San Diego" | as.character(FISCAL_YEAR) == "2003") %>% 
  select(ENTITY_ID, FISCAL_YEAR)

# create the stacked sample
main1 <- data.frame()

for (i in 1:nrow(treated)) {
  treat <- treated[i,]
  if (treat$FISCAL_YEAR < 2011) {
    control <- treated %>% filter(FISCAL_YEAR >= treat$FISCAL_YEAR + 5)
    temp <- rbind(treat, control) %>% 
      mutate(event = i,
             eventyear = treat$FISCAL_YEAR,
             fire = 1*(FISCAL_YEAR == eventyear)) %>%
      select(-FISCAL_YEAR)
    temp <- munipanel %>% left_join(temp, by = "ENTITY_ID")
    temp <- temp %>% filter(FISCAL_YEAR >= eventyear-5,
                            FISCAL_YEAR <= eventyear+4)
    main1 <- rbind(main1, temp)
  }
}

# setting up the regression
main1 <- main1 %>% 
  mutate(eventtime = -1*(fire == 0) + (FISCAL_YEAR - eventyear)*(fire == 1))
main1$eventtime <- relevel(factor(main1$eventtime), ref = 5)
main1 <- main1 %>% 
  mutate(treatDID = 1*(as.character(eventtime) %in% c("0", "1", "2", "3", "4")))

rm(temp, treat, control)


pop <- felm(ihs(POP) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, data = main1)

genrev <- felm(ihs(GENREV_TOT) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, data = main1)

funcrev <- felm(ihs(FUNC_TOT) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, data = main1)

exp <- felm(ihs(TOT_EXP_TOT) ~ treatDID | ENTITY_ID + FISCAL_YEAR | 0 | event, data = main1)

excrev <- felm(identity(GENREV_EXC_OF_GEN_REV/POP) ~ treatDID | ENTITY_ID + FISCAL_YEAR + event | 0 | event, data = main1)

stargazer(pop, genrev, funcrev, exp, excrev, type = "text", out = "table10.txt",
          title = "Robustness checks with alternative exposure measure",
          dep.var.labels = c("Population", "Gen. Rev.", "Func. Rev.", "Total Exp.", "Tot. Exc. Rev."))

rm(main1)



############## Table 11 ##############

main1 <- main %>% mutate(treatDID2 = treatDID*(ENTITY_ID %in% c(1266, 1597, 1183, 1395, 1150)),
                         treatDID = treatDID - treatDID2)

pop <- felm(ihs(POP) ~ treatDID + treatDID2 | ENTITY_ID + FISCAL_YEAR | 0 | event, 
            data = main1)

genrev <- felm(ihs(GENREV_TOT) ~ treatDID + treatDID2 | ENTITY_ID + FISCAL_YEAR | 0 | event, 
               data = main1)

funcrev <- felm(ihs(FUNC_TOT) ~ treatDID + treatDID2 | ENTITY_ID + FISCAL_YEAR | 0 | event, 
                data = main1)

exp <- felm(ihs(TOT_EXP_TOT) ~ treatDID + treatDID2 | ENTITY_ID + FISCAL_YEAR | 0 | event, 
            data = main1)

excrev <- felm(identity(GENREV_EXC_OF_GEN_REV/POP) ~ treatDID + treatDID2 | ENTITY_ID + FISCAL_YEAR + event | 0 | event, 
               data = main1)

stargazer(pop, genrev, funcrev, exp, excrev, type = "text", out = "table11.txt",
          title = "Heterogeneous Effects Based on Wildfire Severity",
          dep.var.labels = c("Population", "Gen. Rev.", "Func. Rev.", "Total Exp.", "Tot. Exc. Rev."),
          covariate.labels = c("Other Fire X Post", "Severe Fire X Post"))

rm(main1)
